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ABSTRACT 



Accurate knowledge of the effect of feedback from galaxy formation on the matter dis- 
tribution is a key requirement for future weak lensing experiments. Recent studies using hy- 
drodynamic simulations have shown that different baryonic feedback scenarios lead to signif- 
icantly different two-point shear statistics. In this paper we extend earlier work to three-point 
shear statistics. We show that, relative to the predictions of dark matter only models, the ampli- 
tude of the signal can be reduced by as much as 30 — 40% on scales of a few arcminutes. As is 
the case for two-point shear tomography, the interpretation of three-point shear statistics with 
dark matter only models is therefore plagued by a strong bias. However, we find that baryonic 
feedback affects two- and three-point shear statistics differently and demonstrate that this can 
be used to assess the fidelity of various feedback models. In particular, upcoming surveys such 
as Euclid can discriminate between different feedback models by measuring both second- and 
third-order shear statistics. Because it will likely remain impossible to predict baryonic feed- 
back with high accuracy from first principles, we argue in favour of phenomenological models 
that can capture the relevant effects of baryonic feedback processes in addition to changes in 
cosmology. We construct such a model by modifying the standard (dark matter only) halo 
model to characterise the generic effects of energetic feedback using a small number of pa- 
rameters. We demonstrate that weak lensing surveys such as Euclid may be able to mitigate 
the effects of baryonic processes, such as outflows driven by feedback from star formation 
and AGN, by marginalising over the feedback parameters. 

Key words: Gravitational lensing:weak, surveys - large-scale structure of the Universe - 
cosmological parameters - Cosmology:theory - galaxies:formation 



1 INTRODUCTION 

Weak lensing by large-scale structure, i.e. cosmic shear, is sensitive 
to the matter distribution without having to rely on assumptions 
about its nature or dynamic state. Thus, it represents a powerful 
tool to study the statistical properties of matter density fluctuations. 
In particular, the measurement of the two-point shear statistics as 
a function of the redshift of the sources, allows one to study the 
evolution of the power spectrum of the matter fluctuations (e.g. 
Hu 2002). In this respect, weak lensing is complementary to pri- 
mary cosmic microwave background (CMB) experiments, which 
by themselves do not constrain the late-time evolution of density 
fluctuations. For a recent review of cosmological applications, see 
Hoekstra & Jain (2008). 

Weak lensing measurements from future surveys such as 
c/j43 ( [Laureijs et al. 2009 ; Amen dola et al. 2012| l, will offer the 
possibility to test the validity of the standard cosmological model 



and, in particular, to discriminate between dark energy and mod- 
ified gravity models. They will have the statistical power to con- 
strain the equation of state of the dark energy at the percent level. 
However, a change of a few percent in the equation of state cor- 
responds to a change of a percent in the amplitude of the matter 
power spectrum. This means that in order to exploit the statistical 
power of those datasets, one needs to predict the power spectrum 
with percent level accuracy up to scales of A: ~ 10 /i/Mpc for 
which the evolution of the fluctuations is in the non-linear regime 
(e.g. Hearin et al. 2012). 

Currently, the interpretation of two-point shear statistics re- 
lies on predictions which are based on N-body simulations of col- 
lisionless particles. To date, this has been sufficient, as the mat- 
ter present in the Universe is mainly coUisionless and statisti- 
cal uncertainties large. However, baryons represent about 17% 
( IKomatsu et al. 20Tll of the total matter content in the Universe 
and the power spectrum of density fluctuations therefore also de- 
pends on their distribution. 
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The distribution of baryons is governed by physical mecha- 
nisms such as radiative cooling, star formation and feedback from 
supemovae and active galactic nuclei (AGN). Those mechanisms 
have been implemented only recently in hydrodynamic simulations 
and the accuracy of the results is still under discussion. Springel 
(2010a, 2010b, 2011) thoroughly discuss the performances of the 
smoothed particle hydrodynamics (SPH), adaptive mesh refine- 
ment (AMR) and moving mesh methods. Along these lines, Vo- 
gelsberger et al. (20 11) demonstrated that different numerical tech- 
niques can lead to different properties on galactic scales. However, 
Schaye et al. (2010) and Scannapieco et al. (2012) suggest that the 
differences between numerical techniques are smaller than the dif- 
ferences due to variations in the physical mechanisms when con- 
sidering the properties of groups. In particular, McCarthy et al. 
(2010), Puchwein et al. (2010), Fabjan et al. (2010), and Dubois 
et al. (2011) have shown that the presence of AGNs significantly 
affects the properties of groups, leading to much improved agree- 
ment with observations. 

Clearly, an accurate description of the matter distribution in 
the Universe cannot be achieved without a correct description of the 
baryonic component. Indeed, van Daalen et al. (201 1) showed that 
baryonic feedback processes modify the power spectrum of matter 
fluctuations (see also Jing et al. 2006; Levine & Gnedin (2006); 
Rudd et al. 2008; Guillet et al. 2010; Casarini et al. 2011) and 
that, depending on the scenario, this modification can affect sur- 
prisingly large scales. In particular, van Daalen et al. (201 1) found 
that in their model with AGN feedback, their only model for which 
predictions for groups are in good agreement with both X-ray and 
optical observations (McCarthy et al. 2010), the matter power spec- 
trum is suppressed by 1 per cent and 10 per cent at ~ 0.3 and 
1 ft/Mpc, respectively. Semboloni et al. (2011a) used those same 
power spectra to derive the cosmic shear signal for different sce- 
narios and showed that ignoring baryonic feedback when interpret- 
ing two-point shear statistics can strongly bias the measurements 
of cosmological parameters. While these studies focused on two- 
point statistics, higher order statistics must also be affected. For 
instance, Yang et al. (2012) recently investigated how the presence 
of baryons affects peak statistics. For their study, they tried to re- 
produce the steepening of the halo density profiles in the central 
regions caused by star formation, by manually increasing halo con- 
centrations of dark matter only simulations. They find a mild in- 
crease of high peaks, but the effect could be different when using 
hydrodynamic simulations. 

In this paper we will investigate the effect of baryonic feed- 
back on three-point shear statistics, which have been advocated in 
the past as an additional cosmological probe. For instance, previous 
studies have shown that the combination of two- and three-point 
shear statistics provides much tighter cosmological constraints 
than the two-point shear statistics alone (Bernardeau et al. 1997 ; 

ITakada & Jain 2003al .Kilbinger & Schneider 2005 ( 

|Berge et al. & 20I0[ IVafaei et al. 20I0I [Kayo et al. 2012^ . The 
benefit of combining three- and two-point shear statistics is even 
larger in the self-calibration regime ( IHuterer et al. 2006t . where 
one exploits the different sensitivity of the two measurements to 
the source redshift distribution in order to alleviate the impact 
of the lack of accurate knowledge of the redshifts of the source 
galaxies. 

A few attempts have been made to measure three-point shear 
statistics JPen et al. 2003llJarvis et al. 20041 1. Those early measure- 
ments used small datasets and shear measurement algorithms 
which were much less accurate than the ones available today. While 
in both studies a signal was detected, the results were affected by 



residual point spread function (PSF) systematics. More recently, 
Semboloni et al. (2011b) used high-quality space-based data, the 
COSMOS-HST dataset, to perform a measurement of three-point 
shear statistics which does not show evidence of residual systemat- 
ics. 

For large datasets, such as the ongoing Kilo-Degree Survey 
(KiDS), and future surveys such as the Dark Energy Survey (DES), 
the Large Synoptic Survey Telescope (LSST) and Euclid, the sta- 
tistical power of three-point shear statistics is comparable to that 
of two-point shear statistics ( IVafaei et al. 20I0l l. Thanks to its dif- 
ferent dependency on the cosmological parameters, three-point 
shear statistics will thus be a complementary test for cosmolog- 
ical models. However, the enhanced sensitivity to the non-linear 
features of the field of density fluctuations makes the interpreta- 
tion of three-point shear statistics challenging. To date, this inter- 
pretation has relied on models based either on perturbation theory 
( Scoccimar ro & Frieman 19991 Scoccimar ro & Couch man 20()TJ 
or on the halo model ( Taka da & Jain 2 003a; Taka da & Jain 2003ab . 
These models have only been tested against dark matter simulations 
and do not agree to better than 10 — 20%. Moreover, the effect of 
baryonic physics on three-point shear statistics, or equivalently, on 
the bispectrum of matter fluctuations is still unknown. 

In this paper we use hydrodynamic simulations from the Over- 
whelmingly Large Simulations (OWLS; Schaye et al. 2010) to 
assess the impact of including baryonic feedback on higher or- 
der statistics. The paper is organised as follows: in Section [2] we 
describe the method used to compute three-point shear statistics; 
in Section [3] we briefly describe the simulations; in Section |4] we 
present the results from the OWLS models. In Section|5]we intro- 
duce a new parametrisation for the effect of baryonic feedback on 
the matter power spectrum which updates the one by Semboloni et 
al. (2011a). In Section[6]we perform a likelihood analysis show- 
ing how the simultaneous measurement of the two- and three-point 
shear statistics can detect the presence of feedback. In this way, 
three-point shear statistics can be used to self-calibrate for the (un- 
known) baryonic feedback. Furthermore, we show how a general 
parametrisation of the effect of feedback, the strength of which is 
only marginally constrained, can be used to reduce the effect of the 
bias on the cosmological parameter constraints. Finally, we discuss 
our results and conclude in Section[7] 



2 CALCULATION OF THREE-POINT SHEAR 
STATISTICS 

In this section we recall briefly the relation between cosmic shear 
and density fluctuation statistics and introduce the notation used in 
this paper. We define the matter density contrast 5 such that the den- 
sity field at a position x is p{x) = p(l + (5(a;)), with p the average 
density. We define the power spectrum P(k) and the bispectrum 
B{ki,k2, fca) as: 



5{kl)S{k2)j = (27vySoikl+k2)P{kl), (I) 

5{ki)S{k2)S(k-a)) = {27vfSuiki + fea + fc3)B(fci, fca, fcg) ,(2) 



where 5 is the Fourier transform of the density contrast 5 and k 
is the Fourier conjugate of x. The fact that ^5(fci)5(fe2)^ differs 
from zero only if fci = — fea and depends only on the modules 
of the wave vector is a consequence of assuming 5 to be a sta- 
tistically homogeneous and isotropic field. This also implies that 
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{^5{k\)5{k2)5{k->,)j is defined only for closed triangles and does 
not depend on their orientation. 

The effect of foreground overdensities on background galaxies 
is a stretching (i.e. a circular source would appear elliptical) and a 
(de)magnification. Those effects are completely described by the 
complex shear 'y = 71 + 172 and the convergence Both are 
related to the second order derivatives of the projected gravitational 
potential (see for example Schneider et al. 1998). In particular: 



K S) = - 



dwg{w)fk{w) , (3) 

a(w) 



where s is a two-dimensional vector such that s = k±fk{'w) with 
k± the projection of k perpendicular to the line of sight. Hq is 
the Hubble constant, Qm is the matter density parameter, c is the 
speed of light, w is the comoving radial coordinate, fkiw) is the 
comoving angular diameter distance at radial distance w, a is the 
cosmological scale factor. 



p^[w )dw — , , (4) 



and pni {■w)d'w describes the distribution of sources as a function of 
the radial coordinate, or equivalently as a function of the redshift, 
in which case it is denoted as pz {z)dz. 

Analogous to Equations (TJ and we define the convergence 
power spectrum and bispectrum: 



(k(si)k(s2)) = (27r)25D(si + S2)P4s) , 



(5) 



{«;(si)k(s2)k(s3)) = {2tt)^5o{si + S2 + S3)-Bk(si, S2, S3) .(6) 

One cannot measure n directly, but one can instead estimate the 
shear 7 by measuring the ellipticities of galaxies. In fact, assum- 
ing that galaxies are randomly oriented, any observed correlation 
between the orientations of galaxies is the result of the deflection 
of light by foreground structures and thus a direct measurement of 
the shear correlations. The power spectra of k and 7 are the same 
and the relation between two-point shear correlations and Pk{s) is 
straightforward. 

The relation between the three-point shear correla- 
tion functions and the bispectrum is more complicated 
dSchneider et al. 2005t . For this reason, the measured three- 
point shear correlation functions are generally combined to obtain 
a new third-order statistic: 

(A4(e)) = j d^eiUe{0i) j d'e2Ue{e2) j d^e^Ueie-i) (?) 



I 



(27r)2 



(27r) 



(k(si)k(s2)k(s3)) 



which is directly related to the projected bispectrum. In fact, by 
choosing 

the relation between three-point correlation functions and 
(Mfp(6l)) is relatively simple l lPen et al. 2003] IJarvis et al. 20041 
[Schneider et al. 2005b . For completeness, we remind the reader 
that: 



— (k(si)k(s2)) . 



(9) 



(27r)2 J (27r) 
Inserting Equations ^ and ^ into Equation ([7} and using the 



definition in Equation one derives an explicit relation between 
the 3D bispectrum of density fluctuations and (Mfp(0)). Thus, to 
measure (i\/fp(6')) for various sets of simulations, we in principle 
need to estimate the bispectrum. Measuring the 3D bispectrum di- 
rectly in each simulation as a function of redshift is very time con- 
suming and we leave this for future work. Here, we make use of 
an approximation derived using perturbation theory ( [Fry 1984] ) and 
adapted to cold dark matter (CDM) cosmologies by Scoccimarro & 
Frieman (1999) and Scoccimarro & Couchman (2001): 

B{ki,k2, k'i) = 2F(fei, fc2)P(fci)P(fc2) + cycl. (10) 

where cycl. indicates a cyclic permutation of the indices, 

5 

F(ki,k2) = -a{n,ki)a{n,k2) (11) 

1 feife2 /fci k2\ . , , , 
+ --r-;— ( T + y jK^^y fci)fe(n, fe) 



+ 



2 fcifc2 V k2 ki 

2 /feifeo\2 

7 



where n is the spectral index of the primordial power spectrum. 
The coefficients 



a(n, k) — 
h(n, k) — 
c{n, k) = 



1 + a8(^)"°-'[0.7Q3W]^/^(g/4)"+«-" 
1 + (q/4)"+»-5 

1 + 0.4(n + 0.3)g"+^ 

1 + qn+3.5 

l + 4.5/[1.5 + (n + 3)''](2g)"+3 



1 + (2g)"+3-5 
have been fitted to N-body simulations with 

4 — 2" 



(12) 



(13) 



(14) 



(15) 



1 + 2"+i ' 

a-s{z) is the standard deviation of density fluctuations linearly 
evolved to redshift z, integrated in a sphere of radius 8 Mpc/h. The 
variable q = k/k„i with 47rfcjJiPiinoar(fcni) = 1 defines a remap- 
ping. Scoccimarro & Couchamnn (2001) showed that this approxi- 
mation is accurate to within 15% up to fc of a few /i/Mpc. 

Semboloni et al. (2011b) compared (Affp(^)) results from 
dark matter only simulations to the above approximation and the 
non-linear halofit dSmith et al. 2003 l l power spectrum. They found 
that (Mfp(6')) is systematically underestimated on small scales, 
although the underestimation may be caused by the fact that the 
halofit approach underestimates the power spectrum at small scales 
(e.g. Hilbert et al. 2009, Semboloni et al. 201 lb) and not necessar- 
ily by the perturbative expansion. Even in the case in which the 
perturbative expansion breaks down and does underestimate the 
value of the three-point cosmic shear statistics, the relative com- 
parison between different feedback scenarios does not suffer from 
this limitation and is expected to be accurate as long as the form of 
F{ki, k2) is not affected significantly by the presence of baryons. 
We argue that this is the case, because the existence of baryonic 
feedback should not affect the transition from the linear to the non- 
linear regime. In future work we will investigate the accuracy of 
the approximations used in this paper by directly measuring the 
bispectrum in the simulations. 



3 SIMULATIONS 

We focus here on the four simulations from the Overwhelm- 
ingly Large Simulation (OWLS) project ( jSchaye et al. 20I0l l 
that were also analysed in Semboloni et al. (2011a): 
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k (h/Mpc) k (hIMpc) k (h/Mpc) 

1.0 0.1 1.0 0.1 1.0 0.1 

1.1 . . • • n 1.1 ' ' ' ' n 102 n ' ' ' ' ^ 




0.5 ^ , , , , J 0.5 ^ , , , , J 0.96 , , , , J 

1 10 1 10 1 10 

9(arcmin) d(arcmin) 9(arcmin) 

Figure 1. The left panel shows the ratio of (A/fp(9)) of the REF (green solid line), DBLIMV1618 (pink soHd line) and AGN (blue solid line) results and 
the DMONLY values as a function of the angular' scale 9 for a survey with a similar depth as Euclid . The dashed Hnes show the same results for a KiDS-like 
survey. The middle and right panels show the same as the left panel but for ^Af|p(6')) and S-i{8), respectively. We indicate on the top the physical scale k 
that contributes most to the value of (A/^p(6)) and (Affp((?)) assuming a source redshift of z = 1. 



DMONLY, REF, DBLIMFV1618 and AGN. All the sim- 
ulations have been realised using the same initial con- 
ditions from a ACDM WMAP3 l |Spergel et al. 2003) 
fiducial cosmology: {Qm,Qb,^A,cTs,n,h} = 
{0.238,0.0418,0.762,0.74,0.951,0.73}. The length of the 
box side is L = 100 Mpc The DMONLY simulation con- 
tains 512^ collisionless particles with mass 4.1 X 10* Mq h-^. 
The hydrodynamic simulations contain the same number of 
collisionless particles plus an equal number of baryonic particles 
with mass 8.7 x 10^ Mq . We describe here the other main 
features of these simulations. 

• DMONLY: a dark matter only simulation, of the kind com- 
monly used to compute the non-linear power spectrum for weak 
lensing studies. It is therefore the reference to which we compare 
the other simulations. 

• REF: this simulation represents a standard scenario assumed 
in cosmological hydrodynamic simulations. It includes most of the 
mechanisms which are thought to be important for the star forma- 
tion history (see Schaye et al. 2010 for a detailed discussion) but 
not AGN feedback. The implementation of radiative cooling, star 
formation, supernova driven winds, stellar evolution and mass loss 
have been described in Wiersma Schaye & Smith (2009), Schaye 
& Dalla Vecchia (2008), Dalla Vecchia & Schaye (2008), and 
Wiersma et al. (2009) respectively. 

• DBLIMFV1618: this simulation has been produced using the 
same mechanisms as REF. The only difference between the two 
simulations is that in this simulation the stellar initial mass func- 
tion (IMF) was modified to produce more massive stars when the 
pressure of the gas is high, i.e. in starburst galaxies and close to 
galactic centres. This is obtained by switching from the Chabrier 
(2003) IMF assumed in the REF model to a Baugh et al. (2005) 
IMF in those regions. The IMF change increases the number of 
supernovae and thus enhances the effect of stellar winds suppress- 
ing the star formation rate (SFR) at small redshifts. However, this 
mechanism alone is not able to reproduce the observed SFR (see 
Schaye et al. 2010). 

• AGN: a hydrodynamic simulation which differs from REF 
only by the inclusion of AGN. The AGN feedback has been mod- 
elled using a modified version of the prescription of Springel et al. 
(2005) as described in Booth & Schaye (2009). In this approach 



the AGN transfers energy to the neighbouring gas, heating it up 
and driving supersonic outflows which are able to displace a large 
quantity of baryons far from the AGN itself. 

In addition, the models WDENS, NOSNJSfOZCOOL and NOZ- 
COOL, which are described below, are used in Section[5]to test the 
generality of our feedback model, but their impact on the cosmic 
shear signal is not evaluated. 

• WDENS: a hydrodynamic simulation which differs from REF 
in that the initial mass loading and velocity of the supernova driven 
wind vary as a function the surrounding gas density so that they 
become more efficient way to transfer kinetic energy even in high 
pressure regions. 

• NOSNJSfOZCOOL: The SN feedback is removed and the ra- 
diative cooling rate is computed assuming primordial abundances. 

• NOZCOOL: Cooling rate assumes primordial abundances. 

We know from van Daalen et al. (2011) that the modifications to 
the power spectrum of the WDENS simulation are quite close to 
the ones of DBLIMFV1618 thus we expect deviations in cosmic 
shear signal which are also similar. The NOSNJSfOZCOOL and 
NOZCOOL scenarios are unrealistic and/or cannot reproduce a re- 
alistic SFR. The power spectra in these two last simulations are not 
modified at large scales (see van Daalen et al. 201 1), thus we expect 
changes in the cosmic shear signal similar to the ones of the REF 
scenario. 

Among the hydrodynamic simulations considered here, the 
AGN is arguably the most realistic, as it is able to reproduce the 
gas density, temperature, entropy, and metallicity profiles inferred 
from X-ray observations, as well as the stellar masses, star forma- 
tion rates, and stellar age distributions inferred from optical obser- 
vations of low -redshift groups of galaxies ( [McCarthy et al. 20lOp . 



4 RESULTS 

In this section we evaluate the effect of different baryonic feed- 
back models on the predicted value of (A/fp(^)). We use the mat- 
ter power spectra P{k,z) measured by van Daalen et al. (2011) 
together with the approximations l|10t-(|15t to estimate the bispec- 
trum. Finally, we integrate the bispectrum over the source distri- 
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bution using Equations ([3}, @ and Q. We chioose thie following 
source distribution: 



zoV{- 



(16) 



where , Pz and zo have been fixed to simulate a KiDS-like survey 
and a Euclid-like survey. Following Vafaei et al. (2010), we use 
the following values: = 0.81, Pz = 3.15, zo = 1.19 for the 
KiDS-like survey and = 0.96, = 1.70, zo = 1.07 for the 
Euclid-like survey, corresponding to median redshifts of 0.8 and 
0.91 respectively. 

We show in the left panel of Figure[T]the ratios of (i\ffp(&)) 
for REF, DBLIMFV1618 and AGN to {Mip{9)) for DMONLY as 
a function of the angular scale 9. Note that the aperture-mass fil- 
ter defined by Equation ([8} is fairly broad and for a given angular 
scale 9 peaks at s = V2/9 (where 9 is given in radians). As a con- 
sequence, aperture mass statistics depend strongly on the density 
field at small scales. As an example, for a typical weak lensing sur- 
vey, the value of (Affp(S)) at 6* ~ 1 arcmin depends on the power 
spectrum P{k) at fc of a few ft/Mpc (see Figure 2 of Semboloni et 
al. 201 la). Similar to what was found by Semboloni et al. (201 la), 
the suppression is in this case about 20% at arcminute scales for the 
AGN scenario, whereas it is much smaller for the REF scenario. 
The DBLIMFV1618 is an intermediate case between the REF and 
AGN cases. 

We show in the middle panel of Figure [T] the ratios for 
^ A'/fp (^) ) . The ratio is close to unity for the REF scenario, whereas 
it has a large-scale dependence in the case of AGN feedback, and 
the amplitude of (A/fp(S)) is suppressed as much as 40% on ar- 
cminute scales. Moreover, the impact is significant on large scales, 
being about 10% at 61 ~ 10 arcmin. 

Finally, the right panel of Figure[T]shows the effect of baryonic 
feedback on the skewness: 



S3{0) 



{Mi,{9)Y 



(17) 



The skewness has been considered a very appealing measurement 
as it scales cx and depends only weakly on erg. In principle, it 
is an effective measure to break the parameter degeneracy typical 
of two-point shear statistics (Bernardeau et al. 1997). However, it 
has recently been shown IVafaei et al. 20lOl l that it is it much less 
effective to constrain cosmology even when combined with two- 
point shear statistics. Indeed, the skewness has a fairly flat profile 
at small scales where the two- and three-point shear statistics show 
features. As a result, these statistics are more sensitive to cosmol- 
ogy. The right panel of Figure [T] shows that the skewness is also 
relatively insensitive to the baryonic feedback scenario (note the 
smaller y-axis range). One of the goals of this paper is to investi- 
gate whether one can detect the effect of feedback by measuring 
weak lensing statistics. For this reason, we focus here on the mea- 
surement of (A/fp(6l)), which is more sensitive to variations in the 
feedback model than the skewness. 



shown by Seljak (2000) (see also Mandelbaum et al. 2005), the halo 
model reproduces the power spectrum into the non-linear regime 
quite well, although some parameters have to be calibrated using 
numerical simulations. The power spectrum is computed as the sum 
of two terms: 



(18) 



The first term describes the correlation of the density fluctuations 
within the same halo. This Poisson term P^ (k) dominates on small 
scales and is given by 



(27r)3 



duf{u)^y[k,Miu)f, 



(19) 



where p is the mean matter density and y[k, M{i')] is the Fourier 
transform of the density profile of a halo with virial mass M{iy) 
which is normalised such that: 



y[k,M] 



/J^^" 4:nr'^drp{r) 



(20) 



where p(r) is the density profile of haloes truncated at the virial 
radius r^ii. The virial radius is defined in practice as the radius 
within which the average density is greater than Svh times the mean 
density of the Universe. We compute S^ir as in Henry (2000) and 
Nakamura & Suto (1997). The peak height i/ of such an overdensity 
is defined as 



Sciz) 



a{M) 



(21) 



where Sc is the linear theory value of a spherical overdensity which 
collapses at a redshift z and a{AI) is the rms of fluctuations in 
spheres that contain a mass M at an initial time, extrapolated to 
z using linear theory. To compute the mass function f{iy) we use 
dSheth & Tormen I999I ): 



i^f{u) = A{1 + v'-^y^l'^ exp(-i.72), 



(22) 



where v' — ai/ with a = 0.707 and p = 0.3. The normalisation 
constant A is determined by imposing / f{i/)dv = 1. Note that the 
function f{iy) is related to the halo mass function dn/dM through 



^dM = -^f{iy)du. 
dM M-'^ ' 



(23) 



The second term, P^^{k), describes the clustering of haloes and 
dominates on large scales. It is given by 



P''''(fc) = PUk)[ / duf{u)h{u)y[k,M{u)] 



(24) 



where Piin(fc) is the linear power spectrum, and the halo bias h{i/) 
is given by JMandelbaum et al. 2005b : 



5 TOWARDS A GENERAL DESCRIPTION OF 
FEEDBACK 

5.1 Halo model and baryonic feedback 

Throughout the paper we use the halo model to predict the matter 
power spectrum analytically. The reason for this choice is that the 
halo model framework allows one to naturally account for the mod- 
ifications to the power spectrum caused by baryonic feedback. As 



+ 2p ^25) 

with a = 0.73, p = 0.15 and v' = av. To predict the weak lensing 
signal one needs to integrate Equation l |19t and l l24t in a range of 
masses that is wide enough. We take 10** < A/vir < IO^^Mq/Ii. 
Larger haloes have a very small probability to form whereas small 
haloes do not contribute significantly to the weak lensing signal. 
The last ingredient that we have not yet specified is the density 
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Figure 2. Left panel: relation between the halo mass A/500 ™d the total baryonic mass contained within rsoo of the halo. The baryonic mass content that 
corresponds to the universal baryonic fraction as a function of A/soo is indicated by the black solid line. The dashed lines show the best fit power laws to the 
corresponding points. The solid lines show the best fit to the points when we impose the constraints that all the scenarios reach the universal baryonic fraction 
for a halo the same given mass A/pivot which is about 2.7 X 10^^ A/q/Zi ( more details about its deiivation are given in the text). Right panel: relation between 
the hot gas and the stellar mass within A/500 . The dashed lines represent the best-fit power laws to the data. 



profile of the matter haloes. Numerical cold dark matter simulations 
have shown that the NFW profile ( INavarro et al. 19931 l 



p{r) oc 



(26) 



where is the scale radius, is a fair description of the radial matter 
distribution for haloes with a wide range in mass. They also indi- 
cate that Ta is not a free parameter, but that it is related to the virial 
mass (albeit with considerable scatter). It is customary to account 
for this correlation by specifying a relation between the concentra- 
tion c = r^ir/rs and the virial mass. When baryonic physics is 
included, the dark matter profiles remain fairly well approximated 
by NFW profiles (see for example Duffy et al. 2010), although the 
mass-concentration relation changes due to the back-reaction of the 
baryons on the dark-matter component. 

We have already shown in Semboloni et al. (201 la) that mod- 
ifying the traditional NFW profile in the halo model to account for 
gas profiles and star fractions provides a good description of the 
power spectrum predicted in each feedback scenario. A somewhat 
different approach has been developed by Zentner et al. (2008) 
who account for the eifect of feedback by changing the mass- 
concentration relations. The advantage of the latter approach is that 
halo profiles are still parametrised with only two parameters al- 
though the mass concentration relation changes as a function of 
redshift. However, if the hot gas and the dark matter have a very 
different distribution this approximation might not be very accu- 
rate. Furthermore, it is more difficult to constrain using observa- 
tions than the model developed by Semboloni et al. (2011a) that 
we employ here. 

We assume that each halo is collapsing with an initial frac- 
tion of baryons which is equal to the universal gas fraction, /gas = 
flb/ Qrn ■ We compare this initial fraction to the fraction of baryons 
(hot gas plus stars) contained in rsoo, defined as the radius of a 
sphere within which the mean density equals 500 times the criti- 
cal density. We select the same haloes as Semboloni et al. (201 la) 
which are also used in McCarthy et al. (2010). These haloes cover 
a mass range 5 X IO^^Mq/H < M500 < lO^Mo//!. This covers 
the regime of halo masses to which cosmic shear is most sensitive 



(see for example Semboloni et al. 201 la). In this range, the relation 
between the total mass of baryons contained within rsoo versus the 
total halo mass M500 is fairly well described by a power-law (see 
left panel of Figure |2]l: 



logio(Mb 



logio(^^oo) + 6b a 



(27) 



A power law can also be used to describe the relation between the 
hot gas (i.e. gas with temperature T > 2.2 x 10^ K ) and the stellar 
mass contained in haloes (see right panel of Figure|2]l: 



loglO (A/star) = fflatar logio(Mgaa) + fas 



(28) 



For each scenario we use these relations to determine the stellar 
and hot gas fractions contained in haloes. We fit these functions to 
the results from McCarthy et al. (2010, 2011). We extrapolate the 
results for lower and higher masses which are needed for the halo 
model. Although this extrapolation at low masses might not be ac- 
curate it will affect the weak lensing predictions only marginally as 
small masses (M500 < W'^^ Mq / h) do not contribute significantly 
to the cosmic shear signal. When we extrapolate to high masses we 
reach a mass for which the baryonic fraction is equal to the univer- 
sal value. We know from the observations that for higher masses the 
gas fraction does not exceed the universal fraction and we therefore 
fix it to this value. 

We describe the spatial distribution of each of the three com- 
ponents in the halo using different profiles. We approximate the 
stellar profile by a point mass. For the dark matter component 
we use an NFW profile. The hot gas is described by a /3-model 
( rCavaliere & Fusco-Femiano 1976) 

r / r \ 2-1 -3/3/2 
Pgas^PO 1+ (29) 

L Varsoo/ J 
where the parameter a describes the relation between rsoo and the 
characteristic scale of the gas profile, and /3 is the characteristic 
slope of the profile. The total mass of the gas within rsoo deter- 
mines the value of po. The values of a and /? are measured in the 
simulations. The gas that has been moved beyond the virial radius 
is distributed uniformly in a region between rvir < r < 2rvir. This 
choice is partially motivated by independent studies that suggest the 
presence of gas in the outskirts of haloes, but it is mostly a practical 
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Figures. The top-left panel shows the amplitude of (Affp(e)) as a function of the angular scale 9 for the DMONLY (black), REF (green), DBLIMFV1618 
(pink) and AGN (blue) simulations compared with the amplitude predicted using our modified halo model (dashed lines). The top-right panel shows the same 
results as the left panel but for (Mfp(0)). The middle panels show explicitly the ratio between the value measured in the simulations and the predictions from 
the modified halo model and the bottom panels show the same ratio once that the common inaccuracy of the halo model has been accounted for (see text). 
We indicate on the top the physical scale k that contributes most to the values of (Affp{9)) and (A/^p(9)) assuming that the sources are at redshift 2 = 1. 
The plots show the good agreement between the cosmic shear signal measured for any of the scenarios and the predictions from the extended halo model 
developed in this paper. In particular, the middle and bottom panels show that this model describes the cosmic shear signal in any of the feedback scenarios 
with an accuracy that is comparable to the accuracy with which the traditional halo model describes the signal in a dark matter only scenario. 



choice allowing us to modify halo profiles without inserting a new 
component. The validity of this assumption should be investigated 
further, but it does not have a major impact on the conclusions of 
this paper. 

The final profile is then the sum of the individual components: 
j/tot[fc,M] =^j/,[fc,M]^ (30) 

2 = 1 

where i identifies the component of the halo profile, which can be 
dark matter, hot gas within the viral radius, stars or the smooth 
gas component beyond the virial radius. The sum of the masses of 
the four components is the virial mass. The profiles yi[k, M] of 
the dark matter and hot gas components are obtained by truncating 
the density profiles at the virial radius as in Equation l l20t . The 
expression j/i [k, M] for the smooth gas component beyond the viral 
radius is obtained by changing the integration range in Equation 

((20) to Tvir < r < 2rvir. 

To predict the DMONLY model, we use the halo model with 
the usual NEW profile and the mass-concentration relations de- 
rived by Duffy et al. (2008). Note that Duffy et al. (2008) used a 
WMAP5 cosmology which is different form the OWLS WMAP3 
cosmology. Since the mass-concentration relation depends on the 
cosmology (on ag in particular, e.g. Eke et al. 2001) this can intro- 
duce a difference between the measured and predicted spectra. In 
the case of baryonic feedback we use the mass-concentration mod- 
ifications derived by Duffy et al. (2010) using the OWLS models 
in the following way: the differences between feedback models and 
the dark matter only scenario as a function of redshift are obtained 
by linear interpolation of the average of the ratios shown in that 



paper at redshifts z = and 2 = 2. For DBLIMFV1618 Duffy 
et al. (2010) do not provide the mass-concentration relations and 
therefore we chose to use the ones that are derived for a strong SN 
feedback model ("ZC_SFB" in Duffy at al. 2010) which is denoted 
as WDENS here. The a and /3 parameters are also fixed to the av- 
erage value measured in the simulations. Eor DBLIMEV1618, for 
which we could not measure these parameters directly, we fix these 
values to the AGN scenario as we found that this provides us with 
a better model than when using the REF profile. 

In Figure [3] we show the signal measured in the simulations 
compared with the predictions obtained using the halo model, both 
for (A/fp(^)) and {Al'ip{0)) for a survey like Euclid. Since we are 
using the wrong mass-concentration relations and we are fixing the 
stellar fraction, the gas fraction and the gas profiles at 2; = 0, we do 
not expect our model to work perfectly. Nevertheless, the model can 
reproduce the general features for all feedback scenarios quite well. 
In particular, our approach reproduces (A/fp(&)) and (A/fp(S)) 
predicted by the simulations to within 10% and 20% respectively, 
for any of the baryonic feedback descriptions. This accuracy is sim- 
ilar to the one with which the dark matter only halo model describes 
the DMONLY results. This is shown more explicitly in the middle- 
left panel where we plot the ratio 

„_ (MfpW°wL^) 
(Mlp^hod) ' 

i.e. the ratio of the values measured in the simulations and the ones 
predicted by the extended halo model. In the right-middle panel, 
we show the same results for (Affp(6l)). 

The lack of accuracy of our model at intermediate scales. 
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about 10% for (Affp(6l)) and 20% for {Ml^{e)) at scales of 
^ ~ 10 arcmin, is mainly caused by the lack of agreement 
between the halo model and the DMONLY results. The lim- 
ited accuracy with which the halo model in its original formu- 
lation predicts dark matter power spectra has been discussed in 
a number of papers (see for example l IHilbert et al. 2009l >) and 
it is beyond the scope of this paper to improve the model or 
to discuss its performance. It is, however, clear that its perfor- 
mance needs to be improved before it can be used to inter- 
pret data from large weak lensing surveys. The "halo-exclusion" 
dTinker et al. 20051 ICacciato et al. 2009t and in particular the per- 
turbative halo model approach ( [Valageas & Nishimichi 20rT| l are 
notable steps in this direction. Note, however, that the discrepancy 
between the DMONLY and halo model prediction is not necessar- 
ily due only to the inaccuracy of the halo model. The DMONLY 
results are affected by cosmic variance (see appendix Al of van 
Daalen et al. 2011). Furthermore, as we have mentioned already, 
the mass-concentration relations we have used are not the most ap- 
propriate for the WMAP3 cosmology assumed in the simulations. 

To avoid any issues originating from the direct comparison be- 
tween DMONLY and halo model predictions, we will follow Sem- 
boloni et al. (2011a) and work only with ratios. The advantages of 
using ratios is that we can suppress the effect of cosmic variance 
and finite size but explore the effect of the baryonic feedback up to 
small k. For this reason we divide for each scenario the ratio R by 
the same ratio obtained for the DMONLY scenario and show these 
results in the bottom panels. In this way, one can see better how 
well the model captures the modifications introduced by our feed- 
back model. The corrections we have applied to the halo model re- 
produce the weak lensing signal from the REF and DBLIMFVI618 
scenarios very well, whereas the effect of baryonic feedback for the 
AGN scenario is slightly underestimated. However, the difference 
is only a few percent for (A/fp(6)) and ten percent in the worst 
case for {M%{e)). 

5.2 A general description of feedback 

The model we have described in the previous subsection has 
enough flexibility to adequately describe any of the OWLS sce- 
narios we have examined. In this model, we have assumed that the 
effect of baryonic feedback can be captured by six parameters: two 
describing the shape of the gas profile (a and /?), two describing 
the total mass of baryons (stars -I- hot gas) contained in a halo of 
mass M500, and two describing the stellar mass as a function of 
the mass of the hot gas. Furthermore, for each simulation we used 
the mass-concentration relations from Duffy et al. (2010). Since 
the real strength of the feedback is not known, one could imagine 
marginalising over these parameters to obtain unbiased cosmolog- 
ical constraints. However, one would like to keep the number of 
nuisance parameters as small as possible, both for computational 
and signal-to-noise purposes. The aim of this section is to check 
whether it is possible to reduce the number of degrees of freedom 
in this model and to see if in doing so we are able to describe the 
impact of feedback generically. 

As a first step, we show below that the fraction of stars and the 
fraction of hot gas as a function of halo mass are correlated. In fact, 
looking at Figure[2l one can see that for any given feedback model 
there is a relation between halo mass, mass of hot gas and mass of 
stars; if the feedback is stronger there is less gas and there are also 
fewer stars. This suggests that for a given set of parameters describ- 
ing the gas fraction there may be a set of parameters which defines 
the stellar mass. Furthermore, the halo mass for which the fraction 



of baryons inside the halo equals the universal fraction is very sim- 
ilar for all scenarios. The new fit performed, fixing the 'pivot' mass 
based on the AGN model (logj^Q (Mpivot) = 14.43), is shown in the 
left panel Figure|2](solid lines) and can be compared with the origi- 
nal fit performed for each scenario individually (dashed lines). The 
fit for the REF scenario is unchanged, while the DBLIMFVI618 is 
slightly tilted, but the quality of the fit is the same and the differ- 
ences are minor. Replacing the expression of A/pivot in Equation 
( |27b we obtain: 

logio(A'/b (Afsoo) -logio(A/pivot)] (32) 

+ loglo(A^pivot) + logio(/gas) 

which assigns to any halo of Afsoo a universal gas fraction if 
flbar = 1, and a lower baryonic mass for a given halo of mass 
Afsoo as the slope Obar increases. As before, for masses larger than 
Afpivot we keep the gas fraction fixed to the universal value. We 
cannot reduce the number of parameters for the stellar fraction, 
which we continue to parametrise with Ostar and 6atar. However, 
Figure |4] shows that there is a linear relation between the slope 
Obar and the parameters bstar and Ostar. That suggests that inde- 
pendently of the mechanisms that remove the gas, there is a tight 
relation between star and gas content provided that a single process 
is responsible for expelling the gas and suppressing the star forma- 
tion. Figure |4] also suggests that this relation can be expressed in 
a rather simple way. In order to check the validity of our relation, 
we also show the results of the best-fit obtained for the WDENS 
simulation. Although the linear relations in Figure|4]have been ob- 
tained using only the simulations REF, AGN and DBLIMFV1618, 
the WDENS best-fit values seem to agree well with this relation. 

As we anticipated, the WDENS simulation has a power 
spectrum modification which is very similar to the one of 
DBLIMFV1618 (see van Daalen et al. 2011). However, the phys- 
ical mechanisms responsible for changing the large scale matter 
distribution are different and it is remarkable that even in this case 
the gas and stellar content seem to follow the same law that we 
found for the other simulations. This suggests that there might be 
some physical reason behind our findings. 

The correlation between the parameters breaks down for the 
NOZCOOL and NOSN_NOZCOOL simulations. Whereas it is true 
that simulations which do not implement metal enrichment and 
metal-line cooling and simulations without any energetic feedback 
from star formation are unrealistic, it is not necessarily true that 
including these processes allows one to obtain simulations with 
a realistic star formation history (the REF simulation is an exam- 
ple). However, in the simulations where metal-line cooling and su- 
pernova feedback are included the conditions for star formation 
are similar and the injection of energy in the gas followed by its 
displacement determines how frequently these conditions are met. 
This common mechanism seems to be responsible for the relations 
we have found. Whether these simple relations are able to describe 
a large class of realistic feedback scenarios is a question that needs 
to be addressed with other sets of simulations using a wider variety 
of feedback implementations and is beyond the scope of this paper. 

If these relations were proven to be generic, they could provide 
us with a recipe to determine the hot-gas and the stellar fractions 
within a halo of a given mass. This would mean that the feedback 
could be described by a single parameter: the fraction of gas left in 
a halo as a function of its mass. The parameters of the /3-model and 
mass-concentration relation characterise the overall distribution of 
matter within the haloes and also change the predicted power spec- 
trum. In principle, they also depend on the specific feedback model, 
but we expect them to be fairly degenerate. For this reason, we keep 
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Figure 4. Solid red triangles (blue solid circles) show the value of Ostar 
(bstar) describing the fraction of stars in each halo (see Equation |28K as 
a function of the total baryonic content expressed by the value of a^ar 
(see Equation |32) . From left to right the values of a^ar correspond to REF, 
DBLIMFV1618, WDENS and AGN. The error bars provide the errors from 
each individual fit. Note that if abar = 1 the baryonic fraction is the same 
as the universal gas fraction of any halo. The dashed lines are the best linear 
fits to the data from REF, DBLIMFV1618 and AGN. These fits can capture 
the general effect of feedback on the stars and hot gas distribution in the 
haloes also for the WDENS simulation. Open diamonds indicate the results 
for the NOZCOOL and NOSN.NOZCOOL scenarios, which do not follow 
the trend displayed by the other models. 



a, 13 as free parameters in our model when we perform a likelihood 
analysis in the next section. The ranges are chosen to include val- 
ues measured in the simulations and are a = [0.0005,0.15] and 
/3 — [0.45, 0.85]. We do fix the mass-concentration parameters to 
those for the AGN simulation. This is just one possible approach 
to marginalise over the lack of knowledge of the overall profile. 
One could also fix the /? profile and change the NFW parameters or 
vary them both within the ranges set by the simulations. Note that 
the NFW parameters do not change much between the various feed- 
back scenarios, while the /3 profile parameters do. The dispersion in 
those parameters is also large and depends on the mass. However, 
since the hot gas fraction is small, the modifications to the overall 
profile have only a minor impact on the power spectrum: to first 
order all the modifications arise because of the removal of the hot 
gas. 



5.3 Cosmology dependence 

Finally, we investigate the cosmology dependence of the feedback 
model adopted in this paper. To do so we compute the ratio of 
power spectra P{k, z) / (k, z) predicted by our halo model 
for different cosmologies. Whereas, the universal gas fraction /gas 
fixes the fraction of gas present in each halo at the moment of 
collapse, Equation i32l determines which fraction of this initial 
mass is still within the halo at low redshift and which fraction is 
ejected beyond the virial radius. For this reason any change in fim 
leads to the same modification of the power spectrum as long as 
/gas ~ const. Note that, since we do not know how the profile 
parameters vary as a function of cosmology, we keep them fixed. 

When we change the normalisation of the matter power spec- 
trum (78, our halo model predicts that the ratio of the power spec- 
tra changes slightly. This differs from van Daalen et al. (2011) 
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Figure 5. Ratio between the dark matter power spectra predicted by our 
halo model for various baryonic scenarios and the dark matter only case, 
for 2 = and different cosmologies. The WMAP3 cosmology results (solid 
lines) are compared to the same results obtained for a WMAP3 cosmology 
but with (Tg = 0.84 (dashed lines). In this case /gas is the same as the 
WMAP3 cosmology and the differences are small. The dot-dashed lines 
show the results for a WMAP3 cosmology but with /gas = 0.120. In this 
case the ratio in the case of DBLIMFV1618 and AGN scenario is signif- 
icantly changed. The solid thin lines represent the ratio measured in the 
simulations (which use the WMAP3 cosmology) to show the accuracy of 
the model used in this paper. 



who have shown that, in the case of the AGN scenario, the rela- 
tive variation to the dark matter only power spectrum is very close 
for the WMAP3 and WMAP7 cosmologies. We remind the reader 
that these cosmologies differ mainly because of the change of erg . 
The discrepancy between our results and the ones in van Daalen 
et al. (2011) can be understood as follows. Changing the value of 
o-g changes the halo mass function. One can see from Equations 
l ll9t and l l24b that the halo mass function acts as a weight applied 
to any profile of mass Al. The ratio between the total power spec- 
tra is given by the ratio of the weighted sums. Since we keep the 
profile parameters unchanged when we change the cosmology, it 
is not surprising that the ratio varies. On the other hand, we know 
that the profile characteristics vary with the cosmology. In partic- 
ular, smaller values of ag lead to the formation of haloes with a 
smaller concentration, at least in DMONLY scenarios. It remains to 
be seen if our model is able to reproduce the result by van Daalen 
et al. (201 1) when the profiles are assigned correctly as a function 
of cosmology, something which we avoid doing here because we 
do not know how the profiles vary as a function of cosmology, es- 
pecially when baryonic feedback is included. 

We show in Figure [5] the ratio of the power spectrum in each 
scenario to the dark matter-only scenario for two different cosmolo- 
gies, and compare these with the WMAP3 results. As one can see 
comparing the solid lines (WMAP3 cosmology) and dashed lines 
(WMAP3 cosmology with different ag value), varying as while 
keeping the universal baryonic fraction constant leads only to small 
variations of the ratio. When the universal gas fraction is changed 
(dot-dashed lines in Figure[5]l the modifications to the power spec- 
trum are also changed. The variation is not large in the case of REF 
scenario as the fraction of gas that leaves the haloes is small, but 
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it is larger in the case of the AGN and DBLIMFV1618 scenarios. 
In particular, when the universal gas fraction decreases, the effect 
of the feedback is reduced as the dark matter fraction in the halo is 
increased. 

Finally, we note that this model is different from the one we 
used in Semboloni et al. (2011a) because we use different scaling 
relations. In particular, we now fit the stellar masses of the halos 
directly, whereas Semboloni et al. (201 la) fitted the luminosity and 
used observed mass-to-light ratios to derive the stellar fraction. Be- 
cause of the changes and the dispersion in the fitted relations, which 
for simplicity were ignored in both papers, one expects the model to 
perform slightly differently, although this does not change the con- 
clusions of Semboloni et al. (201 la). The description adopted here 
is more general than the one adopted in Semboloni et al. (2011a) 
as Equation l |32t provides the general relation needed to modify the 
dark matter power spectrum in the way we have just shown. 



6 LIKELIHOOD ANALYSIS 

In this section we will show that the combined measurement of 
two- and three-point shear statistics on large datasets can unveil the 
existence of strong baryonic feedback. Furthermore, we will show 
that the effect of an unknown feedback scenario can be mitigated, 
provided it can be modelled as described in the previous section. 
To do that, we perform a likelihood analysis in the [S7m,o"8, wo] 
space. The posterior probability is given by: 

P(7r|d) oc r(d|7r)P(7r) , (33) 

where: 

£(d|7r) = exp[-l/2(d - m(7r))*C"^(d - m(7r))] (34) 

and we have indicated with tt, d and m the parameter space, the 
data and the model vectors, respectively. To compute the covari- 
ance matrix C, we use the set of projected density maps described 
in Vafaei et al. (2010) which have been realised with the same 
WMAP3 cosmology. Each convergence map consists of 1024^ pix- 
els covering a total area of 12.84 deg^. We have 60 lines-of-sight 
and along each line-of-sight the projected density field is sampled 
in 40 steps between < 2 < 3. For any line-of-sight, we com- 
bine the 40 projected convergence maps using the redshift distri- 
bution of each survey we are interested in. Finally, we add shot 
noise to the 60 convergence maps assuming a galaxy density of 
10 galaxies/arcmin^ for the KiDS-like and 30 galaxies/arcmin^ 
for the Euclid-\\ks survey. The intrinsic ellipticity dispersion is set 
to 0.2 per component. We use an area A — 15000 deg^ for the 
Euclid-like survey and ^ = 1500 deg^ for KiDS. We use 12 an- 
gular scales in the range 0.5 to 20 arcmin. The maximum angular 
scale is then smaller than the side of each map and we can rescale 
the covariance matrix by multiplying by the ratio of the area of the 
maps to the area of the survey we want to simulate. In fact, for a 
given survey depth and galaxy density, the overall statistical error 
scales oc A'^^^, assuming that boundaries effects can be neglected 
and that the survey is composed of contiguous observations (see for 
example Schneider et al. 2002). Strictly speaking this is only true in 
the Gaussian regime, and since we apply this simple rescaling our 
covariance matrix is understimated (see Sato et al. 2009 and Kayo 
et al. 2012). Although it is important for the future to either simu- 
late very large volumes or to use analytic formula as suggested by 
Kayo et al. (2012), this is not a major issue for this paper. 

To compute the predicted cosmic shear signal, we use the lin- 
ear power spectrum of density fluctuations derived with the transfer 



function from Eisenstein & Hu (1998). The non-linear evolution 
is derived using our halo model. The data vector d(^) is derived 
computing (Mfp(e')) and (Mip{9)) for each OWLS scenario. To 
alleviate effects arising from the lack of agreement between halo 
model (dark matter-only) predictions and the DMONLY measure- 
ment, we follow Semboloni et al. (2011a) and we recalibrate the 
data vector d(^): 

where d{9)'~'^^^ is the data vector measured in the baryonic 
OWLS simulations and d(6')°^°'^^'^ is the data vector measured 
for the dark matter only simulation. This ensures that when the in- 
put vector is the one measured in the DMONLY scenario, the data 
vector we use in the likelihood analysis is the one obtained using 
d(6i)''°'*, the halo model with the same cosmology as OWLS. 

We show in the left panel of Figure [6] the la error contours 
in the f2m-o"8 parameter space inferred from the measurement of 
(J\ffp(^)) (solid contours) and (Affp(^)) (dashed contours) for the 
Euciid-like survey. We have used a flat prior on wq — [—1.2, —0.8] 
which we assume can be inferred from non weak-lensing experi- 
ments. The left panel of Figure |6] shows that the interpretation of 
both statistics from the Euclid-like survey, using a dark matter only 
model, provides a probability distribution of the posterior which is 
increasingly biased when the feedback becomes stronger. Further- 
more, there is a discrepancy between the posterior probability dis- 
tributions for (A/fp(6')) and {Mip{6)). The consequence of bary- 
onic feedback is to displace the posteriors in a direction more or less 
perpendicular to the — Qm degeneracy direction, suggesting that 
the simultaneous measurement of two- and three-point statistics 
can be used to detect strong baryonic feedback. The stronger the 
feedback, the more gas is removed from the inner part of the haloes, 
leading to larger differences between the posteriors for (Mfp(&)) 
and for (A/^p(&)). In particular, one can see in the right panel of 
Figure[6]that strong feedback models, such as AGN, could already 
be tested, albeit at the la level by measuring two- and three-point 
shear statistics with the KiDS survey. Note that the difference be- 
tween this figure and Figure 6 of Semboloni et al. (201 la) is caused 
by the fact that we use here statistics and scales which are more sen- 
sitive to the baryonic feedback. In general dark matter only models 
are a bad description of baryonic feedback but the amplitude of the 
bias in any parameter depends on the scales used to perform cosmic 
shear analysis and on the priors. However, there is another factor 
which might cause small differences between posteriors generated 
in this paper and Semboloni et al. (2011a): in that paper non-linear 
predictions were computed using the fitting formula of Smith et 
al. (2003) whereas here we used the halo model. Both approaches 
suffer of limited accuracy and the cosmological parameters depen- 
dence is not the same. This issue has been thoroughly discussed 
by Hearin et al. (2012) who also quantified the difference of the 
likelihood contours using different methods. 

It is interesting to investigate whether the interpretation 
of {M'ip{e)) from the CFHTLS-Wide (Hoekstra et al. 2006J 
IFu et al. 20081 [Heymans et al. 20T2) dataset is already sensitive to 
baryonic feedback. To mimic the CFHTLS-Wide dataset we as- 
sume the same source distribution as Euclid , a survey area A = 
170 deg^, a galaxy density n = 15 gal/arcmin^ and an intrinsic 
ellipticity dispersion a^ ~ 0.3 per component. Under these as- 
sumptions, we find that (i\ffp(S)) may be biased at the la level 
in the case of strong feedback, such as the AGN scenario, and if all 
cosmological parameters but Q,,n and erg are fixed. However, the 
statistical uncertainty for (Mfp(^)) is too large for the two like- 
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Figure 6. Left panel: likelihood analysis results for (Mfp(6»)) (solid lines) and (Mfp(e)) (dashed lines) for the DMONLY (black), REF (green), 
DBLIMFV1618 (pink) and AGN (blue) scenarios. The Icr contours have been obtained assuming Euclid-like fiducial survey. The square indicates the WMAP3 
cosmology which should be the maximum of the probability distribution in case of unbiased results. Right panel: same as the left panel but for a KiDS-like 
survey. For the KiDS-like survey we show only deviations for DMONLY and AGN scenarios. The models do not include baryonic feedback and have been 
computed using the halo model power spectrum and rescaled as described in the text. Any shift relative to the DMONLY case indicates the presence of bias 
due to baryonic effects. 



lihood posteriors to be in tension. Finally, we note that (7l/ap(S)) 
does not capture all the information contained in the three-point 
shear statistic jSchneider et al. 20051 1. For example, by measuring 
(Map(6'i)Map(6'2)Map(03)) One is able to capture better the sta- 
tistical power of the three-point shear statistics. We did not include 
such measurements here as we have only a few lines-of-sight to 
estimate the covariance matrix. We expect that this would increase 
the discrepancy between two- and three-point results by reducing 
the errors on the posteriors of as and Qrn from three-point shear 
statistics. This might also be the case if one uses tomography. 



6.1 Calibration of baryonic feedback 

We have shown in the previous section that the occurrence of bary- 
onic feedback can be unveiled by a simple likelihood analysis of 
two- and three-point statistics measured on small scales where the 
effect is expected to be largest. In the case of strong feedback a dark 
matter only interpretation will lead to a discrepancy of the posterior 
distributions inferred from the two types of statistics. 

We have also shown in Section [5] that the effects of energetic 
feedback in each of the simulations can be described by a generic 
model that allows us to parametrise the strength of the baryonic 
feedback irrespective of the underlying physics that is responsible 
for removing the gas from the central parts of the haloes. More- 
over, we have argued that there is a relation between the gas con- 
tent and the stellar content in the simulations. This allowed us to 
parametrise the strength of the feedback by a unique parameter Obar 
which changes from one scenario to another. 



We will now use this parametrisation to study whether it is 
able to describe a general feedback model in such a way that the fi- 
nal cosmological analysis for a generic unknown feedback scenario 
will be unbiased. We vary the parameter abar in Equation i32\ be- 
tween the extreme values measured in the OWLS simulations. For 
a given Obar, we will assign the baryonic content and the stellar 
content to the haloes according to the best-fit relations shown in 
Figure |4] 

As we already discussed, the mass concentration and the gas 
profile parameters depend also on the feedback scenario (and on 
the cosmology). Since these parameters should be degenerate, we 
fix mass-concentration relations to the one of the AGN simulation 
and keep a and /3 as free parameters. In this way we have three 
nuisance parameter plus three cosmological parameters. 

As we previously mentioned, it has it has been shown by van 
Daalen et al. (20 11) that the change from a WMAP3 to a WMAP7 
cosmology does not change the suppression of power caused by 
the baryonic feedback. We will assume here that this is also the 
case for other changes in the cosmological parameters, although we 
have shown that our halo model contradicts this. We have argued 
that this can be explained by the fact that the dark matter profiles 
we use when we change cosmology are fixed to the ones obtained 
for the WMAP5 cosmology plus feedback modifications derived 
within the WMAP3 cosmology. This is an important limitation of 
the current model but it could be overcome when more simulations 
are available. To circumvent this problem we compute the modifi- 
cations for a WMAP3 cosmology and we assume that these are the 
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same for any other cosmology, by computing: 

P{k, zr = P{k, zf^ X g^'^S^^^ (36) 

to obtain P{k, z)™, the estimate of the power spectrum for a given 
cosmology and feedback model. Note that in this way a, (3 and abar 
become nuisance parameters that span the possible range of modifi- 
cations of the power spectrum generated by baryonic feedback, but 
do not necessarily correspond to a physical quantity in a cosmology 
other than WMAP3. We include among the possible scenarios the 
WMAP7 AGN for which we also derive the posterior probability. 

Note that given the small number of simulations we refrain 
from combining the likelihood analysis of two and tree-point shear 
statistics. Instead we would like to show that our model is good 
enough to reduce the bias on the cosmological parameters for both 
statistics. In the left panel of Figure |7] we show the results ob- 
tained for the REF, DBLIMFVI6I8, AGN scenarios and the case 
of a WMAP7 cosmology with AGN feedback. The cosmological 
constraints are now biased by la in the worst case scenario. The 
posteriors are in good agreement, also for the WMAP7 cosmol- 
ogy. In the middle panel we show the posterior for the parameter 
fibar which differs from one scenario to another. As was already 
shown in Figure[5] the modelling is not yet perfect. Therefore, we 
expected to have some residual bias. Note, however, that the cos- 
mological constraints are generally derived using measurements up 
to much larger scales than the ones we used here. The maximum 
angular scale we used in our analysis is 6 = 20 arcmin but current 
surveys already allow measurements on angular scales larger than 
a few degrees and the bias is therefore expected to be smaller in 
practice. 

In the right panel of Figure [7] we compare the posterior prob- 
ability for the NOZCOOL and NOSN_NOZCOOL models that do 
not follow the same empirical law we have applied to establish the 
fraction of gas and stars in our model (see Section [Sjl. The poste- 
riors indicate only a small bias, even though our model is wrong. 
The reason for this is that these models implement a weak feed- 
back (i.e. the fraction of gas expelled outside the core of the haloes 
is small) which changes the power spectrum mostly at small scales. 
Because of that, the uncorrected bias is small, similar to the one 
obtained for the REF scenario. Marginalising over the feedback pa- 
rameters, even using the wrong model, reduces this bias to within 
the la level for NOZCOOL, whereas the bias remains larger for 
the NOSNJMOZCOOL scenario. 



7 CONCLUSIONS 

Measurements of cosmic shear from upcoming weak lensing sur- 
veys have the potential to provide unique constraints on cosmolog- 
ical parameters and the theory of gravity. To interpret these mea- 
surements, models are required that can predict the signal with 
great accuracy. Van Daalen et al. (201 1) have used simulations from 
the OWLS project (Schaye et al. 2010) to show that the baryonic 
physics of galaxy formation, in particular the strong feedback pro- 
cesses that are required to reproduce X-ray and optical observa- 
tions of groups of galaxies, modifies the distribution of matter on 
such large scales that models will somehow have to account for 
these rather uncertain effects. Indeed, in Semboloni et al. (2011a) 
we demonstrated that cosmological constraints based on two-point 
shear statistics will be biased if they are derived within a dark mat- 
ter only framework. 

Here we extended the analysis of Semboloni et al. (2011a) 



to three-point statistics. For a suite of OWLS models, we com- 
puted the three-point shear statistic (Mfp(0)) using an approxi- 
mate model based on perturbation theory (Scoccimarro & Couch- 
man 2001) that takes the non-linear power spectra tabulated by 
van Daalen et al. (2011) as input. We found that scenarios with- 
out strong feedback, which therefore suffer from the overcooling 
that has plagued most models in the literature, predict values of 
{M'tp{9)) that agree with the dark matter only simulation to within 
a few per cent. However, simulations that include either strong stel- 
lar feedback (because of a top-heavy stellar initial mass function in 
starbursts) or AGN feedback, predict strong deviations from the 
dark matter only case on scales 6 < 10'. In the case of AGN 
feedback, which is the only scenario that reproduces observations 
of groups of galaxies (McCarthy et al. 2010, 2011), the value of 
(Mfp(&)) is suppressed by up to 30 — 40% on scales of a few ar- 
cmin. For comparison, the suppression measured for (A/fp(S)) is 
about 20% on the same angular scales. Hence, three-point statis- 
tics are even more affected by baryonic feedback than two-point 
statistics. 

Because of the enhanced sensitivity to non-linear structure and 
the different cosmology dependence, ignoring baryonic processes 
in the interpretation of three-point shear statistics results in a bias 
that differs from the bias affecting two-point shear statistics. As a 
consequence, the posterior probability distributions for cosmolog- 
ical parameters are different for the two types of measurements. 
Thus, the combination of two- and three-point shear statistics en- 
ables tests of galaxy formation models and allows one to check 
whether cosmological constraints inferred from one type of mea- 
surement may be biased due to baryon physics. We find that this 
tension may already be detectable at the ~ la level using upcom- 
ing KiDS data and that it will be highly significant for the inter- 
pretation of Euclid data. These cosmic shear measurements can be 
used along with more direct observations of the gas distribution in 
groups and clusters, including the outskirts, based on X-ray and 
Sunyaev-Zeldovich data, as well as galaxy-galaxy lensing observa- 
tions. This is important, since our lack of knowledge of the physical 
mechanisms governing the overall distribution of baryons makes it 
impossible to construct simulations that can robustly predict the 
density field from first principles. 

We presented a modified version of the halo model of Sem- 
boloni et al. (2011a) and used it to predict the non-linear matter 
power spectrum. The model captures the effects of baryonic feed- 
back by introducing six parameters which characterise the depen- 
dence on halo mass of the hot gas and the stellar mass fractions as 
well as the shape of the gas profiles. All these parameters can in 
principle be measured from observations, but here we determined 
their values by fitting to the results of the simulations (Duffy et al. 
2010; McCarthy et al. 2010, 2011; Semboloni et al. 2011a). For 
each scenario explored here, the model is able to reproduce the 
cosmic shear signal with an accuracy similar to that with which the 
standard halo model reproduces the dark matter only power spectra. 

As the accuracy with which the standard halo model can re- 
produce the dark matter only predictions is insufficient for future 
surveys, and because solving this issue is not the aim of our study, 
we corrected the predictions of our model for the discrepancy be- 
tween the halo model and the OWLS dark matter only simulation. 
This correction also mitigates any issues related to cosmic variance 
due to the finite size of the simulation volume. 

We found that although different scenarios predict different 
halo gas and stellar fractions, the predictions of all simulations that 
include both metal-line cooling and at least some energetic feed- 
back follow the same relations between the amounts of hot gas and 
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Figure 7. Left panel: probability distribution for fim and erg marginalised over too and the feedback parameters a, fi and Obar for a survey similar to Euclid . 
Solid lines represent the 1 cr contours inferred from ^A/^p(9)), while the dotted Hues represent the same contours but infen'ed from ^A/^p(6)). The squares 
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stars within a halo of a given mass and its baryon fraction. We used 
these relations, which are common to all physically plausible sce- 
narios that we investigated, to reduce the number of parameters 
needed to describe the effects of feedback. 

We used this to perform likelihood analyses on simulated data, 
using a version of our halo model that contains three nuisance pa- 
rameters to describe the effect of feedback on the matter distribu- 
tion. Using this model, and marginalising over the three baryon- 
related parameters as well as over the equation of state of the dark 
energy, we found that the posterior probabilities for the cosmolog- 
ical parameters Q.m and ug based on measurements of (i\ffp(6)) 
and (Mfp(^)) on scales between 0.5 and 20 arcmin are no longer 
discrepant, even for a survey as large as Euclid. With this approach 
the bias on the cosmological parameters became smaller than the 
1(7 confidence regions while the overall statistical power was de- 
graded only mildly. Since cosmological constraints are generally 
obtained from data including larger scales, we expect the bias to 
be even less significant, and the same should be true for the loss of 
information caused by the presence of three extra nuisance param- 
eters. 

One hypothesis underlying the use of our modified halo 
model, is that the feedback modifications are independent of the 
cosmology. This hypothesis is supported by the comparison be- 
tween the WMAP3 and WMAP7 cosmologies reported by van 
Daalen et al. (201 1), but it needs to be tested further using a wider 
range of cosmologies and ideally also using simulations that in- 
clude modifications of gravity. The method also assumes that our 
model and the assumed relations between the baryonic mass frac- 
tion and the mass fraction of stars and hot gas are sufficiently flex- 
ible to describe the general modifications of the dark matter only 
predictions expected in realistic scenarios for baryonic feedback. 
This assumption also needs further testing using hydrodynamical 
simulations with a wider range of feedback implementations. 

Although the performance of the model presented in this pa- 
per is still limited, and many of the underlying hypothesis remain 
to be tested, it represents a significant improvement over methods 
that rely solely on dark matter only predictions. In the future the 



assumptions made in this paper will be investigated further and the 
model will be better calibrated using simulations and available ob- 
servations. 
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